Pade approximants for improved finite-temperature spectral functions in the numerical 

renormalization group 
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We introduce an improved approach for obtaining smooth finite-temperature spectral functions of quantum 
impurity models using the numerical renormalization group (NRG) technique. It is based on calculating first the 
Green's function on the imaginary-frequency axis, followed by an analytic continuation to the real-frequency 
axis using Fade approximants. The arbitrariness in choosing a suitable kernel in the conventional broadening 
approach is thereby removed and, furthermore, we find that the Fade method is able to resolve fine details in 
spectral functions with less artifacts on the scale of a; T. We discuss the convergence properties with respect 
to the NRG calculation parameters (discretization, truncation cutoff) and the number of Matsubara points taken 
into account in the analytic continuation. We test the technique on the the single-impurity Anderson model and 
the Hubbard model (within the dynamical mean-field theory). For the Anderson impurity model, we discuss the 
shape of the Kondo resonance and its temperature dependence. For the Hubbard model, we discuss the inner 
structure of the Hubbard bands in metallic and insulating solutions at half-filling, as well as in the doped Mott 
insulator Based on these test cases we conclude that the Fade approximant approach provides more reliable 
results for spectral functions at low-frequency scales of a; < T and that it is capable of resolving shaip spectral 
features also at high frequencies. It outperforms broadening in most respects. 

PACS numbers: 71.27.+a, 75.20.Hr 



I. INTRODUCTION 

Nonmagnetic host metal doped with magnetic impurities 
exhibits low-temperature anomalies in its thermodynamic and 
transport properties known as the Kondo effect^— . The resis- 
tance of such materials has a minimum at finite temperature 
(the Kondo temperature, Tk^- An unyielding mystery for a 
number of years, this "Kondo problem" of the resistance min- 
imum was eventually theoretically explained by J. Kondo by 
showing that the perturbation theory in higher order applied 
to the s-d exchange Hamiltonian includes logarithmically di- 
vergent terms*.. The effective exchange interaction strength J 
is renormalized to larger values as the temperature is reduced: 
magnetic impurities are weakly interacting at high tempera- 
tures, but become strong electron scatterers on the scale of 
Tfc^"—. The system behaves as a renormalized local Fermi 
liquid (FL) at low temperatures'^'^. Accurate results for the 
temperature dependence of the anomalies in the thermody- 
namic properties were obtained by K. Wilson using a novel 
technique, the numerical renormalization group (NRG)^— . 
It consists in reducing the Hamiltonian for the point-like im- 
purity to its one-dimensional form, discretizing the continuum 
of conduction-band electrons in a logarithmic way around the 
Fermi level, rewriting the Hamiltonian in the form of a tight- 
binding chain with exponentially decreasing hopping matrix 
elements, and numerically diagonalizing the resulting Hamil- 
tonian in an iterative way. The NRG results for the thermo- 
dynamics of the Kondo model were later confirmed by the 
Bethe Ansatz approach, which provides an analytic solution 
to the probleroi^— . 

While thermodynamic properties are simple to compute us- 
ing the NR G'"''^''^''^ , dynamical properties (spectral func- 
tions and various dynamical susceptibility functions) and 
transport properties (conductance, thermopower) require sig- 



nificantly more effort22i22. There are two main difficulties. 
The first concerns the actual calculation of the raw spectral 
functions. A number of algorithmic improvements over the 
years have increased the reliability of the method^—, also 
at finite temperatures--. Since the NRG calculations are per- 
formed for tight-binding chains of finite length, the raw spec- 
tral functions are represented as sets of weighted delta peaks, 
and the second difficulty consists in obtaining a smooth con- 
tinuous representation of the spectra. At zero temperature, one 
usually performs spectral broadening using a log-Gaussian 
kernel which is well adapted to the logarithmic discretiza- 
tion gridJi. With high-quality raw results and using the z- 
averaging trick, it is possible to obtain highly accurate final 
results with few overbroadening effect s-^^'-^'^'-^^ . At finite tem- 
peratures, a simple log-Gaussian broadening kernel is not ap- 
propriate for a; < T (in units of h = ks = 1). Instead, on low 
energy scales one should switch to a Gaussian or Lorentzian 
kernel^. Even when smooth cross-over functions are used to 
glue the broadening kernel for w > T with that for uj < T— , 
the resulting spectral functions exhibit artifacts in the cross- 
over frequency region which cannot be fully eliminated. It is 
thus desirable to devise new approaches for obtaining more 
accurate continuous representations of spectral functions. 

Impurity problems have applications ranging from mag- 
netically doped materials and heavy-fermion compounds^ii^, 
electron transport in nanostructures (quantum dotsi^"— , car- 
bon nanotubes^', molecules and single-atom transistors^—, 
molecular magnet a^^-^^) , to spectra of single atoms probed us- 
ing a scanning tunneling microscope (STM)^"— and dissipa- 
tive two-state systems^, but their importance has also vastly 
increased in recent decades because they play a central role in 
the dynamical mean-field theory (DMFT) for strongly corre- 
lated electron materials^—. In the DMFT, a lattice problem 
of correlated electrons (such as Hubbard model, periodical 
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Anderson model, or Kondo lattice model) is mapped onto an 
effective single-impurity problem subject to self-consistency 
conditions. The approach is based on the observation that in 
the limit of infinite dimensions or infinite lattice connectiv- 
ity, the self-energy function S(k, oj) becomes purely localS, 
and depends only on frequency u, not on momentum k. In 
this limit the method is exact, but it is also a good approxi- 
mation for three-dimensional and some two-dimensional sys- 
tems where spatial correlations are less important. Many tech- 
niques have been used in the past to solve the effective im- 
purity models, which is the most computationally demand- 
ing part of DMFT calculations. Presently, the most popular 
methods are exact diagonalization, quantum Monte Carlo, es- 
pecially the continuous-time algorithm (CT-QMC)^"— , and 
NRG. The QMC is a stochastic simulation which is numer- 
ically exact, but time consuming. The results for spectral 
functions are accumulated in bins defined on the Matsubara 
imaginary-frequency axis and to obtain the spectra on the real- 
frequency axis one has to perform an analytic continuation ^i, 
for instance using the maximum entropy method^ or with 
Pade approximants-iaZl. This procedure is not without peril 
and high quality numerical results are required to obtain reli- 
able spectra. The NRG, on the other hand, can provide the 
results directly on the real-frequency axis and is thus well 
adapted to study the very low temperature properties and vari- 
ous details in the spectra. At finite temperatures, the most vex- 
ing technical difficulty of the NRG is the broadening problem 
discussed previously. Since the DMFT consists of iteratively 
solving the effective impurity problem, the spectral function 
artifacts at lj ^ T tend to amplify. This is particularly both- 
ersome for calculations of transport properties^— where the 
contribution to the integrals comes mostly from the frequency 
window uj e [-5r : 5r]. Another issue in the DMFT(NRG) 
approach is the calculation of the self-energy as the ratio of 
two Green's functions 
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((K,gint];4)). 



(1) 



where i^int is the interaction part of the Hamiltonian, while 
is the impurity-orbital annihilation operator The main prob- 
lem is that the causality requirement ImE(aj + iS) < can 
be (slightly) violated at low temperatures around the Fermi 
level. Ironically, the difficulties occur at low temperatures and 
on low-frequency scale, i.e., in the regime which is commonly 
believed to the the forte of the NRG, while it performs rather 
better than expected on intermediate and high temperature and 
frequency scales. 

In this work we explore an alternative approach to spec- 
tral function calculations at finite temperatures. We propose 
to calculate the Green's functions on the discrete set of the 
imaginary Matsubara frequencies and then (when necessary) 
compute the Green's function for real-frequency arguments by 
analytic continuation using the Pade approximants. Switching 
to the imaginary axis in the NRG may seem like a step in the 
wrong direction, but it is motivated by the following observa- 
tions: i) The NRG spectra are not affected by large stochastic 
errors like the Monte Carlo data. The fine details are therefore 
not lost and can be resolved by analytic continuation, ii) One 



directly obtains both real and imaginary parts of the Green's 
function on the imaginary axis, thus Kramers-Kronig trans- 
formations are no longer necessary, iii) In the DMFT(NRG) 
loops, most steps of the calculation can be performed in the 
Matsubara space; an analytic continuation is only necessary to 
compute the hybridization function on the real axis which is, 
however, immediately integrated over to obtain the discretiza- 
tion coefficients, iv) There is less arbitrariness compared to 
the broadening approach. 

The paper is structured as follows. In Sec. |II]we describe 
the basic steps in the DMFT calculations, the NRG technique, 
and our implementation of the Pade approximation. In Sec. Hill 
we test the method on single-impurity problems: we com- 
pute the temperature dependence of the spectral function of 
the resonant-level model and the single-impurity Anderson 
model, and compare various approaches. In Sec. |IV]we ex- 
tend this calculation with a self-consistency loop and study 
the Hubbard model at finite temperatures. We conclude by 
discussing the relative merits of the different techniques. 



II. METHOD 

A. Dynamical mean-field tlieory 

The DMFT is based on the observation that in the limit of 
infinite dimensions k) — > which implies the pos- 

sibility of exactly mapping the lattice problem onto a single 
impurity problem subject to self-consistency conditions'^. 
Let us consider the Hubbard model^"— 

i?Hubbard = '^i'^i^- M)Ckcr'^k(T + U n^^Ui^, (2) 
kcr i 

describing a lattice of sites indexed by i which can be oc- 
cupied by electrons with spin a =t ^nd a =|. Here ek is 
the non-interacting band dispersion, /i is the chemical po- 
tential, and U the Hubbard repulsion. Furthermore, = 
l/\/]V e*'' '''Ckcr, and ni„ — cJ^Cio- is the local occu- 
pancy. 

In the DMFT, the information about the lattice structure is 
fully captured by the non-interacting density of states (DOS) 



Po(e) 



(3) 



Using the Hilbert transform, one can compute the correspond- 
ing Green's function Go{z) in the full complex plane: 



Go(z) 



po(e)de 



(4) 



For common lattices (such as Bethe lattice, 2D, 3D, and 
infinite-dimensional cubic lattice, as well as for flat band) 
there exist analytic closed-form expressions for Gq{z). 

The Hubbard model maps upon the single-impurity Ander- 
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son model (SIAM)22r.22i8i-H: 



ka 



(5) 



where the operators fka- annihilate an electron in the contin- 
uum bath and rf^ on the impurity site, while n = + with 
n„ = d^dij- The (complex) hybridization function A, defined 

as 



A(z) 



Z- £k 



(6) 



fuUy characterizes the coupling between the impurity and the 
continuuml^, i.e., there can be different (but physically fully 
equivalent) descriptions in terms of the bath energies and 
hopping constants Vk- 

The quantity of main interest is the momentum-resolved 
Green's function 

GkaW-((cka;cL». (7) 
and the corresponding spectral function 

^kff (w) = - -ImGka {uJ + i5). (8) 

TT 

They provide information about the electron dispersion in the 
presence of interactions and can be used to determine the ther- 
modynamic and transport properties of the system within the 
linear response theor y^"^i^^'^° . The output from the impurity 
solver is the local impurity self-energy Scr(z), which in the 
DMFT is taken to be equal to the lattice self-energy of the 
correlated electron problem, thus 



Gk.(2) 



1 



(9) 



+ y« - fk - S^(z) 

The local (k-averaged) lattice Green's function is then 

ly I 

Po(e)de (10) 



[z + ^ - T.a{z)] - e 

The DMFT self-consistency condition relates the local lattice 
Green's function and the hybridization function through 



A.(z) = z + /i- (Gi;J_^(z) + I].(z)) 



(11) 



The hybridization function is then used as the input to the im- 
purity solver in the next step of the DMFT iteration. The cal- 
culation proceeds until two consecutive solutions for the local 
spectral function differ by no more than some chosen conver- 
gence criterium. The approach to self-consistency can be sig- 
nificantly accelerated using Broyden mixing, which can also 
be very efficiently used to control the chemical potential ji in 
fixed occupancy calculations^. 



B. Numerical renormalization group 

Wilson's NRG is a non-perturbative numerical renor- 
malization group approach applied to quantum impurity 
problemsiS. The cornerstone of the method is the logarith- 
mic discretization of the conduction ban d'"-^^'-^^'^^ . The in- 
finite number of the continuum degrees of freedom is re- 
duced to a finite number, making the numerical computation 
tractable. The discretization is chosen to be logarithmic be- 
cause in the Kondo problem the excitations from each energy 
scale contribute equally to the renormalization of the effec- 
tive exchange coupling^. The band is divided into slices of 
exponentially decreasing width. 



/„ = [-A-"\-A-(™+i)]D, 



(12) 



for holes and electrons, respectively, with to > 0. Here D is 
the half-bandwidth, while A > 1 is known as the discretiza- 
tion parameter. A complete set of wave-functions is then con- 
structed in each interval the first chosen so that it couples 
to the impurity, while other Fourier components are localized 
away from it; only the first wave-function is retained, while 
all the others are dropped from consideration'". This approx- 
imation becomes exact in the A 1 limit, but the accumu- 
lated experience with the method indicates that it is reliable 
even for very large A despite the seeming crudeness of the 
approximation^. The problem is then transformed to a tridi- 
agonal basis using the Lanczos algorithm, the result being 
a tight-binding Hamiltonian which is known as the hopping 
Hamiltonian or the Wilson chain. A number of improvements 
have been devised over the years. On one hand, it has been 
shown that it is advantageous to perform calculations for sev- 
eral interpenetrating discretization meshes and then average 
the results2i^. Such z-averaging (or twist-averaging) allows 
more accurate calculations at larger values of A since the dis- 
cretization artifacts tend to largely cancel out for aptly chosen 
meshes. On the other hand, modified discretization schemes 
further reduce some systematic errors of the original Wilson's 
approach22i2^i22. These improvements are crucial if the goal is 
to obtain accurate results. 

The NRG calculations in this work have been performed 
with the "NRG Ljubljana" package^, which consists of two 
parts. The first is a Mathematica program which initializes the 
problem by performing the exact diagonalization of the ini- 
tial Hamiltonian in a chosen symmetry-adapted basis, and by 
transforming the operators of interest to the eigenbasis. The 
calculations are performed using a computer algebra system, 
and the input to the program are expressions in the familiar 
second quantization notation24. All quantities are stored in the 
form of irreducible matrix elements, and the Wigner-Eckart 
theorem is used to take into account the symmetry properties. 
The second part of the package is the C++ code which per- 
forms the iterative diagonalization. This consists of adding 
the Wilson chain sites one by one, each time constructing the 
Hamiltonian matrices in all invariant subspaces, diagonaliz- 
ing them, and transforming all the necessary operator ma- 
trices. The full description at step N of the iteration corre- 
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sponds to the effective behavior of the system on the temper- 
ature/energy scale a; TV ~ A^^/^. Since the Fock space grows 
exponentially along the chain, only a small part of the com- 
puted eigenstates are retained after each step. A convenient 
way to perform this truncation is to keep the states up to some 
suitable energy cutoff i?cutoff defined in terms of wjv- Alter- 
natively, at most A'kocp states are retained. The systematic 
error introduced by truncation is small due to the energy-sc ale 
separation propert y of the quantum impurity models: the 
matrix elements between excitations of vastly different energy 
scales are very small. The dynamical quantities are computed 
using the spectral decomposition (Lehmann representation) 
for the eigenstates of the full Hamiltonian. The original ap- 
proach was based on the observation by Sakai et al.— that as 
one proceeds from one step to the next, the lowest few eigen- 
states split due to the interaction with the added shell states, 
while the intermediate lower levels do not show any essential 
change. The intermediate states thus form a good approxima- 
tion of the eigenstates of the Hamiltonian in the infinite-chain 
limit. For problems where the high-energy spectral features 
depend on the low-energy behavior of the system, the spectral 
function has to be computed taking into account the reduced 
density matrix obtained from the density matrix of the low- 
energy fixed-poin t-^ It has been shown that a complete 
basis for the full Fock space of the problem can be defined by 
judiciously using the information from the discarded part of 
the NRG eigenstates''^''''*. This method does not suffer from 
over-counting of excitations and it fulfills the normalization 
sum rule-'^ At finite temperatures, this scheme can be im- 
proved by including the contributions to the density matrix 
from all NRG shells- . This full-density-matrix (FDM) ap- 
proach is currently the most reliable method for computing 
the finite-temperature spectral functions using the NRG. In 
the DMFT applications, one is particularly interested in the 
self-energy S, which may be computed as the ratio^ 



G^zY 



(13) 



where 



= (([d,,iJi,np];4».- (14) 

This is known as the self-energy trick (or S-trick)^, and it has 
recently also been implemented in QMC calculations^ffi. 

Since the hopping Hamiltonian is finite, the computed raw 
spectral functions are represented as a sum of delta peaks: 



at LOj) is smoothed into 



It? \io/ijJj 
52 



(16) 



i.e., a Gaussian function on the logarithmic scale, where b is 
the broadening parameter, chosen depending on the value of 
the discretization parameter A and the number of interleaved 
discretization meshes N^- Peaks sharper than the width of 
the broadening kernel will appear broader than they truly are. 
Typically, the value of b is chosen to be 0.6 or less, but with 
large Nz and small A it can be much reduced, to the point of 
largely eliminating the NRG overbroadening problems (at the 
cost of much longer computation time) ^^'^^ . At finite temper- 
atures, the following broadening kernel has been proposed^: 

K{uj,LUj) ~ L{Lu,LUj)h{u!j) + G{uj,ujj)[l — h{ujj)], (17) 

where 

2' 



G{u},u}j) 

h{u}j) 



1 



exp 



ln|a;/wj| 



■ 7 



1, 

cxp 



■ exp 



OJQ 



log l^j/c^o 



(18) 



, \ < UJQ 



(15) 



Here a is the broadening parameter for the log-Gaussian part, 
equivalent to b in the kernel in Eq. ( fTSI l. 7 = a/4, while luq 
is the cut-off where the log-Gaussian goes smoothly into the 
Gaussian part; typically ujq is chosen to be of the order of 
the temperature T. This broadening approach leads to sizable 
artifacts on the scale of uiq. We find that in practice it is better 
to use slightly modified kernel: 

K{uJ,UJj) = L{LU,LUj)h{Lu) +G{uj,LUj)[l - h{i0)], (19) 

which differs only in the argument of the cross-over func- 
tion h. This breaks the normalization condition, but produces 
smoother spectra and the normalization is reestablished us- 
ing the self-energy trick. This procedure still leads to a slight 
bump around to = wo, which can be further smoothed out 
by averaging over several choices of luq. The artifacts can, 
however, never be completely eliminated. In the following 
subsection we therefore discuss an entirely different approach 
to obtaining a continuous spectral function at finite tempera- 
tures. 



Alternatively, one may use fine-grained binning of these delta 
peaks into very narrow intervals on a logarithmic grid (for in- 
stance 1000 bins per frequency decade). To obtain a mean- 
ingful continuous function, these peaks need to be broadened. 
The original approach to obtaining a smooth curve was by 
Gaussian broadening followed by separate spline interpola- 
tion of the results in odd and even steps, and averaging of the 
two curves21. A better approach is broadening by the log- 
Gaussian distribution function"*": each data point (delta peak 



C. Pade approximation 

The Green's function Ga{z) is related to its corresponding 
spectral function by 



lmG„{u + i6) = -7tA^{uj) 



(20) 



and similarly for ^^(z). Instead of performing the broaden- 
ing, the complex functions Go-(z) and Fa-{z) (both analytic in 
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the upper complex half-plane, cf. Titchmarsh's theorem) are 
evaluated at the Matsubara frequencies 

ioJn = i{2n + 1)ttT (21) 

via the Hilbert transform: 

G.(,,„) = r ^W-Md.,, (22) 

and similarly for Fa-{iujn). Note that there is no need to per- 
form the Kramer-Kronig transformation to calculate the real 
parts of G'(j(w) and Fcr{uj). The equation (l22T l is actually just 
a finite sum 



Ga{iuJn) = -. ^ 

1!,),-. — / 



(23) 



If the z-averaing is used, it is performed at this point: 



n— 1 J ^-J 



(24) 



The self-energy is calculated on the imaginary axis in the same 
way as on the real axis via Eq. ( fT3] l which holds in the whole 
upper complex half-plane. From the self energy in the Mat- 
subara frequencies ^^(iwn), one calculates the local Green's 
function Gioc,cr(j<^Ti) using Eq. ( fTOt . To complete the DMFT 
loop, one needs to calculate the hybridization function on the 
real axis. We first calculate 

A^(icj„) = icj„ + A* - (Gi;J^^(iw„) - E^(iw„)^ . (25) 

This function is analytic in the upper complex half plane so 
we can perform an analytic continuation and evaluate A just 
above the real axis to obtain the (real) hybridization function 



= -Im[A,(c^ + iJ)] 



(26) 



An analytic continuation is also needed to determine the spec- 
tral function 

A,{^) = --\^ [Gioc.a (cj + . (27) 

TT 

The generalized mathematical problem of analytic continu- 
ation can be stated as follows: find an analytic function /(z) 
in the upper complex half plane that coincides with calculated 
values on the discrete set of points 



(28) 



where (zj , /j ) are the known point- value pairs of the function. 
The function /(z) must also obey the asymptotic behavior of 
the Green's function [or Acr(z)], namely 



(29) 



There exist two main numerical techniques for the analytic 
continuation: the maximum entropy methodlSi (MEM) and 



Pade approximation^-. MEM is essentially an improved fit to 
available data, taking into account known physical properties 
of the Green's/spectral function, such as sum rules, possible 
symmetries and higher moments of the distribution. In this 
article, we focus on Pade approximation; the comparison of 
MEM and Pade in the context of NRG are yet to be explored. 

The Pade approximation method is based on the assumption 
that /(z) is a rational function 



/(^) 



PO H VVrZ^ 



qo + qiz + • • • + QrZ^ + z 



r+1 ' 



(30) 



where -pj and qj are the unknown complex coefficients to be 
determined. Inserting the Pade approximant ( |30t into equa- 
tions ( |28] l for each point (zj, /j) generates a linear system of 
equations. Defining a vector of unknowns 



X= bo,Pl,---,Pr,'?0,gi,---,'?r], 



a right-hand-side vector 



and a matrix 



( \ zq zl ■■■ fo /ozo 

1 zi zl ■■■ fl /iZl 



V 



the solution for the coefficients is 



= A^b. 



(31) 



(32) 



(33) 



(34) 



We have assumed that the number of points {zj , fj) is 2r. Al- 
ternatively, one can use recursive relations of the continued 
fractions representation of the Pade approximant to evaluate it 
at specific points^. 

Fitting a rational function is a numerically ill-defined prob- 
lem. Let us define ^ as the ratio between the biggest and the 
smallest element in the matrix A. If Zj are the Matsubara 
frequencies, the ratio is approximately 



e= [(4r+l)7rT] 



±r 



(35) 



where we take the minus sign in the power if the base is 
smaller than 1 . In order to invert the matrix A, the numeric 
precision of 2 log2 C binary digits is needed^l, i.e., much more 
than the standard 53-bit mantissa of 64-bit double-precision 
floating point numbers. For this reason, in our implementation 
of the Pade method we use the GNU Multiple Precision Arith- 
metic Library (GMP) for arbitrary precision floating-point nu- 
merics. To solve the linear system of equations, we perform 
Gaussian elimination without pivoting. The process can be 
performed in 0{r^) multiplications. Each multiplication takes 
0{b log 6) CPU cycles, where b is the number of mantissa bits 
in the floating point used (using fast Fourier transform multi- 
plication). 

When the coefficients of the Pade approximant are known 
to high precision, one can approximately calculate the value 
of the Green's function in any point of the upper complex 
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half plane using Homer's polynomial evaluation scheme. The 
most interesting are the values just above the real axis to com- 
pute the hybridization function (w) via Eq. ( l26l l or the spec- 
tral function Acr{i^) via Eq. (|27] |. The Pade approximant can 
be, however, also used for extrapolation when an insufficient 
number of Matsubara frequencies have been computed. It pro- 
vides a good fit to tails even when the asymptotic behavior is 
not yet reached and simple asymptotic 1/z fit on tails does not 
work. 

In most calculations we use Nm — 2r — 350 Matsubara 
frequencies and internal matrix inversion precision of 1024 
mantissa floating bits. In most cases, taking more points does 
not improve the solution; in fact, sometimes taking less points 
helps avoid some issues. The precision of input points plays 
a big role in how good the Pade approximant is, as already 
discussed in Ref. |73j. 

We remark that one can also determine the coefficients of 
the rational function of order r using the data from more than 
Njn = 2r points by least-squares fitting. 



in. RESULTS FOR A SINGLE-IMPURITY PROBLEM: 
ANDERSON MODEL 

We first consider a single impurity in a metal host described 
by the SIAM, now written as 



k(7 



+ 



H.c 



(36) 



where e is the impurity level, while the chemical potential is 
fixed to zero in this section, ^ ~ Q. The spectral function can 
in general be expressed as 



A(w) = --Im, 

TT V + " 



1 



Y.{uj) - A(w) 



(37) 



When the continuum is modelled as a flat band, i.e., a band 
with a constant density of states p in the interval cj e [—1? : 
D] which couples to the impurity with a constant hopping 
Vk = so that the hybridization strength 



(38) 



is a constant, the complex hybridization function A(z) is 
given as 



A(^) 



TT \z/D+l 



or, on the real-frequency axis for —D<lu<D, 

1, fl-uj/D 



A(a; + iS) 



i H In 

TT 



l+uj/D 



(39) 



(40) 



In the limit C/ = 0, the SIAM is equal to the non-interacting 
resonant-level model whose spectral function is given exactly 



by Eq. ( l37b with = 0. For small F, the spectral function 
is well approximated by a Lorentzian 



F/tt 



(w - e)2 + F2 



(41) 



For larger F, the hybridization self-energy effects become siz- 
able, thus the curve shape deviates from a pure Lorentzian 
form, the peak is shifted to 



F fl-ua/D 
cjQ = e in ' 



(42) 



and, in addition, near the band edges lo = ±D the spectral 
function shows additional features: local minima at 



±1?^ 1 



2F 



(43) 



and local maxima at uj^ax given as the solutions of the equa- 
tion 



7r(a;n 



Fin 



1 - UJmax/D 
1 + WmaxZ-D 



0, 



(44) 



located near the band edges. For e = and moderately large 
F, they are approximately equal to 



= ± 



ttD 

1 - exp ( — — 



(45) 



i.e., the maxima are located exponentially close to the band 
edges. Beyond Wmax the spectral function goes to zero at 
u! = ±D in a continuous way, but with a diverging slope. 
Similar behavior is expected in all models where the imagi- 
nary part of the hybridization function drops to zero discon- 
tinuously, so that its real part (Hilbert transform) has a log- 
arithmic singularity. This is the case, for example, for a 2D 
cubic lattice DOS with 



2FD 

A(z) = 'J^KiDyz'), 



(46) 



where K is the elliptic function. The weight of the band-edge 
features is exponentially small for small F, but may be ap- 
preciable for moderate F of the order of the bandwidth, or 
if the main spectral peak is close to the band edge (for in- 
stance, a few times F). Such features typically are not visible 
in the NRG results because of the overbroadening effects, al- 
though in a careful calculation with small broadening width 
and for many values of the twist parameter z, they may be 
observed^. In the other case, when ImA goes to zero con- 
tinuously, the spectral function decreases monotonously near 
the band edges. We note that in non-interacting models, the 
spectral function does not depend on the temperature. 

The NRG calculations in this section are performed with 
the discretization parameter A = 2, with the z-averaging over 
Nz = 16 values of the twist parameter using a modified dis- 
cretization scheme which eliminates the band-edge artifacts 
of the conventional approach^^^. Spectral functions are com- 
puted using the FDM NRG scheme and the S-trick. The trun- 



cation cutoff is set at Er 



[toff 



IOljn unless otherwise noted. 



Flat band with constant F is used in all cases. 
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A. Non-interacting case, U = 




uj/D 



Figure 1 : (Color online) (a) Spectral function of the non-interacting 
resonant-level model obtained with broadening (a — 0.2) and with 
the Pade approximant. The Pade approximant overlaps nearly per- 
fectly with the analytic result, (b) Spectral functions calculated at 
different temperatures T using the Pade approximant (i.e., using the 
sets of Matsubara points iun = i{2n + l)nT with the raw results 
of NRG calculations performed at chosen temperatures T). The an- 
alytic result does not depend on the temperature. Deviations in the 
numerical results occur for T > F. 

Numerical calculation of a spectral function using the NRG 
is a non-trivial test of the method even in the absence of inter- 
actions. Furthermore, even though the spectral function of a 
non-interacting problem does not depend on the temperature, 
numerical approaches generally have more or less severe dif- 
ficulties reproducing this simple fact. In Fig. [T]we therefore 
first analyze the NRG results for the non-interacting resonant- 
level model obtained by the conventional broadening tech- 
nique with a relatively small broadening parameter a = 0.2 
and by the proposed Pade approximant method, and compare 
them with the exact result. The S-trick is obviously not used 
here, because it would trivially produce the exact result. We 
consider a peak centered at w = e = 0, i.e., at the Fermi 
level, where the NRG is said to have very good spectral reso- 
lution. In panel (a) we clearly observe the advantages of the 
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Figure 2: (Color online) (a) Imaginary part of the Green's function 
for the non-interacting resonant-level model (same parameters as in 
Fig-[D evaluated at the Matsubara frequencies for different temper- 
atures. The analytic solution is shown with solid black line, (b) 
Relative errors of raw Matsubara spectral functions (dots) and the 
conesponding Pade approximants (solid lines). The dots quantify 
the intrinsic NRG errors, while the lines measure the full systematic 
error (intrinsic + analytic continuation errors). 



Pade procedure at finite temperatures (here T = F): the Pade 
approximant produces a spectral function which practically 
overlaps with the analytic formula, unlike the broadening pro- 
cedure which produces a severely distorted spectral peak with 
missing spectral weight at low frequencies (the spectral peak 
tails are, however, well reproduced in both approaches; fur- 
thermore, we note that in the T <C F limit, the broadening 
reproduces the exact results to a very good approximation 
on all frequency scales, thus the problems become manifest 
only at finite temperatures). In panel (b) we compare Pade 
approximant results obtained at different temperatures. The 
best agreement with the exact results is obtained for T < F, 
which is expected since for a temperature scale comparable to 
the characteristic physical scales of the problem, the finite set 
of the Matsubara point provides a well matched sampling. For 
T ^ F, the agreement remains good but requires a sufficient 
number of the Matsubara points N,n (this issue is discussed in 
the following). Furthermore, for T = 10~^ we observe some 
artifacts on the flanks of the curve (these are also discussed 
later on). For T ^ F, we find that the spectral height is un- 
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derestimated also in the Pade approximant approach, but less 
severely than in the broadening scheme. This actually indi- 
cates a limitation of the NRG method itself, not of the Pade 
approximation scheme. Even at high T the Green's function at 
the corresponding Matsubara frequencies should contain the 
necessary information to reconstruct the Lorentzian peak, but 
the output from the NRG itself already has sizable systematic 
errors. We explore this question in more detail by plotting 
the raw Green's function on the Matsubara axis in panel (a) 
of Fig.|2l as well as the relative error in panel (b). The exact 
expression on the imaginary axis is (for e = and flat band): 



G(^y) - - 




We find that the systematic error of NRG is below 0.5%. 



— Broadening 




Figure 3: (Color online) Spectral function of the non-interacting 
resonant-level model with large hybridization T/D = 0.3, where the 
behavior at the band edges uj = —D and uj = D may be compared. 

We study the behavior near the band edges using the 
resonant-level model with a very large hybridization T / D = 
0.3, see Fig. [3] The band-edge peaks are completely washed 
out in the broadening approach and the spectral function has 
long tails in the region where it should be equal to zero. The 
Pade approach resolves the band edges, although it is unable 
to correctly resolve the shape and the width of the near-band- 
edge resonances. Since the Pade approximant is a rational 
function, it cannot describe a precipitous drop to zero at band 
edges. The Pade approach also better describes the behavior 
at the Fermi level, which is underestimated when broadening 
is used. 

We now consider the asymmetric case where the spectral 
peak is located at e 7^ 0, i.e., away from the Fermi level. This 
is the situation where the broadening procedure has notorious 
difficulties. The width of the log-Gaussian kernel is namely 
proportional to the frequency, 

Fbr ~ Oi'-^ ~ ae, (48) 

where a is the broadening parameter If the intrinsic width 
of the spectral feature considered is much less than Fbi , sig- 
nificant overbroadening will occur In the example presented 
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Figure 4: (Color online) (a) Spectral function of the non-interacting 
resonant-level model with a peak centered far away from the Fermi 
level, e 3> r. Standard approach produces severely overbroadened 
result since the kernel width for a = 0.2 is Fbr ~ ae = 10"^ The 
Pade method is much more successful, although the peak is slightly 
shifted away from its coixect position, (b) Comparison of the spectral 
functions computed with the Pade approximant for a range of temper- 
atures. Similar degree of agreement is produced for all T considered. 



in Fig. m we take a rather extreme case of T/D ~ 10 ^, 
while T^,/D w 10~2 for a ^ 0.2 and e = 0.05. The 
plot in panel (a) indeed shows that the broadening procedure 
produces a severely overbroadened peak, which is nearly a 
flat line on the scale of the figure. The Pade approximant, 
on the other hand, does not have difficulties reproducing nar- 
row Lorentzian peaks, as expected. We find, however, that the 
peak is somewhat displaced from the correct position, and is 
slightly wider In panel (b) we show that the spectral functions 
depend on the temperature in a non-systematic way, although 
the result is acceptable for all temperatures considered. The 
deviations from the exact results may again be related to the 
raw output from NRG, rather than to the Pade procedure. 

The spectral function computed using Pade approximation 
is not guaranteed to fulfill the normalization sum rule. For 
peak centered at the Fermi level, the normalization deviates 
from 1 to at most a few per mil. For a peak away from the 
Fermi level, the deviations can be sizable, up to several per 
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cent, especially for larger T. In principle, one could reformu- 
late the Pade approximation as a least-squares-fitting proce- 
dure and implement the normalization sum rule (and perhaps 
additional higher-moment sum rules) as a constraint on the 
parameters. This is an idea worth pursuing. 



B. Interacting case, U > 
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Figure 5: (Color online) Spectral functions for the interacting single- 
impurity Anderson model (SIAM) obtained with broadening (a — 
0.2) and by Fade approximation, with and without the E-trick. Two 
sets of parameters with the same T/U ratio are used, thus the effec- 
tive Kondo exchange coupling constant J is the same. For U /D — 1, 
the Hubbard peaks are located close to the band edges, and additional 
spectral features can be resolved at the edge. For U /D = 0.5 the 
Fade results with and without the E-trick overlap nearly perfectly, 
while for U/D — 1 there are some small differences. 

We now turn to the interacting case with finite U. For 
U > ttT and near half-filling, the system is in the Kondo 
regima^M^iiS^. As the temperature is reduced, on the temper- 
ature scale of U the charge fluctuations on the impurity site 
freeze out and the impurity may be described solely in terms 
of its spin degrees of freedo m^- '^i'"-^ ; using the Schrieffer- 
Wolff transformatior^iS^ the model maps onto the Kondo im- 
purity model 



-ffKondo — ^ £fe/k£r/ko- + • S, 



(49) 



where the Kondo exchange coupling constant J is given by 

2r/7r 2r/7r 



pj 



(50) 



S is the impurity spin operator, while s is the conduction-band 
spin density at the impurity position. At the Kondo tempera- 
ture, given byJ^ 



Tk ''^ C/cxp 



1 

"pJ 



(51) 



the spin degree of freedom is screened by the conduction-band 
electrons (the Kondo effect). The two characteristic scales are 
also reflected in the spectral features: the charge fluctuations 
are associated with spectral peaks at ui = e and lu = e + U 
with half-width at half-maximum of 2r, while the screening 
leads to a Kondo resonance^ pinned to the Fermi level, with 
a spectral width of the order of the Kondo temperature. The 
zero-temperature spectral function at the Fermi level is deter- 
mined by the Friedel sum rule 



,105 



^(0) = 



sin S, 



q.p. 



ttT 



(52) 



where S^.p. is the quasiparticle scattering phase shift in the 
local FL theory of the SIAM. If the system is particle-hole 



(p-h) symmetric, i.e., if 5 



e + 



U/2 



0, then Sq,p, = 



tt/2, irrespective of the U/T ratio, thus there is a constraint 
7rryl(0) = 1. In the deep Kondo regime, U/nT > 1, the 
phase shift is close to 7r/2 even if the p-h symmetry is slightly 
broken. Unless F ^ D, we also expect some features near the 
band edges, just like in the non-interacting case, as discussed 
above. At finite temperatures, the Kondo resonance is washed 
out starting at T Tj^. In this work, we use the definition 
of Tk based on transport properties, i.e., Tk is the temper- 
ature where the conductance through a nanodevice described 
by the SIAM is reduced to one half the conductance quantum. 
The conductance and the definition of Tk are discussed more 
thoroughly in Sec. IIII CI 

In Fig.[5]we plot the spectral functions in the p-h symmetric 
case for two different choices of U and F, but with the same 
r/U ratio. The curves obtained by broadening without the 
S-trick are rather typical of standard NRG results: the Kondo 
resonance is well captured, but the Hubbard peaks are signif- 
icantly over-broadened. For U/D = 1, the spectral density 
remains finite even far outside the conduction band due to the 
long tails of the broadening kernel. Furthermore, no features 
near the band edges are detected. 

The curves obtained by the Pade approach are significantly 
improved even without the S-trick. They show all the ex- 
pected features: the Kondo resonance, the Hubbard peaks, and 
(for U /D = 1) the near-band-edge resonances, followed by a 
fast decay to zero outside the band. 

When the E-trick is used, the improvement is dramatic 
in the broadening approach, since the over-broadening is 
strongly reduced. For the U/D = 1 case, the behavior near 
band edges is also improved. In the Pade approach, however, 
the use of the E-trick has very little effect. We thus conclude 
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that when the Ti-trick is not used, the Fade approximant ap- 
proach vastly outperforms broadening, and that the Fade ap- 
proach significantly reduces the need for using the self-energy 
trick. As discussed previously, the S-trick is still very useful 
to restore the normalization of the spectral function to 1 . 
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Figure 6: (Color online) Spectral function of the single-impurity An- 
derson model for different numbers of kept states, A'kcop , in the NRG 
procedure. The artifacts (for instance the dip at a;/D = 0.5 for 
A^'kcop = 500) are truncation-cutoff dependent, thus they stem from 
the raw NRG results. 

However, when using the Pade approximant we also ob- 
serve some anomalies. Because they appear at different lo- 
cations for different NRG truncation parameters, see Fig. 
we argue that they are a direct manifestation of the systematic 
NRG errors that are reproduced too accurately by the Fade 
approximation. The standard broadening approach hides such 
anomalies, but they are, in fact, present in the spectra and 
may be revealed in "high-resolution" calculations with small 
broadening-kemel width. The Fade approximation, however, 
appears to over-emphasize them. Improved results may be ob- 
tained by averaging over different truncation cutoffs^. 

In the broadening approach, there is significant arbitrariness 
in the choice of the shape of the kernel, its width, and the ap- 
proach to handle the w < T and w > T parts in distinct ways. 
In the analytic continuation approach, there is in principle no 
arbitrariness, since the information about the Green's function 
on the set of Matsubara points fully and uniquely determines 
the Green's function in the whole complex half -plane, in par- 
ticular on the real-axis. The way the continuation is performed 
is, clearly, non-unique, but once the Fade approximation ap- 
proach is chosen, the only adjustable parameter is the number 
of the Matsubara points taken into account in the fitting pro- 
cedure. In QMC calculations, it is practice to use sufficiently 
many Matsubara frequencies to reach the asymptotic 1/z be- 
havior on the imaginary axis which corresponds to choosing 
3> U . (Note that the number of points is temperature 
dependent.) In NRG calculation, however, we observe that 
we can recover the tails using less frequencies. The central 
peak is well reproduced using only N,n ~ 50 frequencies and 
Njn ~ 350 is enough to obtain essentially fully converged 
spectral function, as shown in Fig. |7ja). In Fig. [TJb) we plot 
the relative differences between the Fade approximants for 
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Figure 7: (Color online) (a) Spectral functions of the single-impurity 
Anderson model obtained with Fade approximant using the first A^„i 
Matsubara frequencies, (b) Relative difference between the spectral 
functions computed with N,-n = 1000 Matsubara frequencies com- 
pared to those for lower A''^. 

different numbers of Matsubara points. In the Kondo peak 
region the differences are below 0.01, thus the relative error 
contributed by the analytic continuation itself may be esti- 
mated to be well below one percent. From this one may draw 
the conclusion that the main source of error at low frequen- 
cies in this approach is the systematic error of the NRG itself, 
not the analytic continuation procedure. At higher frequen- 
cies the errors are slightly larger; in particular, at the artifacts 
discussed above the error may locally be up to 10%. 

We find that the shift away from the real axis 5 in the ex- 
pression for the spectral function A{ljj) = —{l/iT)lmG{uj + 
iS) does not play an important role, since typically there are no 
poles of the Fade rational function on the real axis (in which 
case it is, in fact, safe to put S ~ 0). For large S, artificial 
broadening is introduced, see Fig.|8] 

Of particular interest is the low-energy part of the spectral 
function, i.e., the Kondo resonance, because it is directly re- 
sponsible for the conductance anomalies observed in magnet- 
ically doped metals and in semiconductor quantum dots. At 
non-zero temperatures, the shape of the resonance is difficult 
to reliably establish due to the broadening problems. Different 
version of the broadening procedure yield improved results 
in some cases, but worse in others. For example, an "opti- 
mized" scheme which produces less artifacts when there is a 
low-frequency spectral peak may lead to more pronounced ar- 
tifacts in the case of a low-frequency spectral dip (and vice 
versa). There is thus no universal choice. The comparison in 
Fig.|9] indicates that the Fade approximant approach does not 
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Figure 8: Spectral functions of the single-impurity Anderson model 
obtained with Pade approach using different shifts from the real axis. 



create bumps near uj ~ ±T, unlike the broadening. In addi- 
tion, it appears that broadening leads to a significant underes- 
timation of the height of the spectral peak at its center Using 
the E-trick in the Pade approach fixes the normalization prob- 
lem. Figure |9] illustrates one of the main results of this work: 
the Pade approximant approach allows to determine the finite- 
temperature spectral functions with less artifacts on the scale 
of w - T. 




Figure 9: (Color online) Close-up view of the top of the Kondo res- 
onance at a temperature greater than Tk- Note the deviation of the 
peak shape from parabolic behavior when broadening is used. 

A function, analytic in the upper complex half-plane, is 
fully determined by its values on a countably infinite set of 
points. When fitting the Pade approximant, one commonly 
takes the values of the function at the Matsubara points ia;„ 
for a finite number N,n of consecutive values of the index n, 
i.e., n = 1,2, ... , N„i. This is, however, not necessary, nor 
always optimal. In particular, when the temperature T is very 
low, the Matsubara points are very dense and the first N,n 
points do not necessarily contain enough information about 




Figure 10: (Color online) (a) Spectral functions of the SIAM at very 
low temperature obtained using Pade approximations which use non- 
consecutive Matsubara frequencies as input: ~ i{2n + l)7rsr. 
The parameter s quantifies the step length. The inset (b) shows 
a close-up on the Kondo resonance, (c) Green's function on the 
imaginary-frequency axis, indicating the Matsubara points taken into 
account in the Pade approximant construction. The black line is the 
asymptotic 1/ z form. 



the high frequency scales, thus it is not possible to reconstruct 
the full spectral function. In these circumstances, we find it 
convenient to increase the spacing to s, so that we use the 
imaginary-frequency points 



i{2n+l)sTTT. 



(53) 



This corresponds to using the Matsubara points for a higher 
effective temperature sT, even though the raw spectral data 
were computed at the actual physical temperature T. The re- 
sults of this procedure are illustrated in Fig. [TO]for the case 
when the temperature is much below the lowest intrinsic tem- 
perature scale of the problem (here the Kondo temperature 
Tk = 4.6 10-3, thus T/Tk = 2 lO"*). For the standard 
choice of consecutive Matsubara points, s = 1, we find results 
which are clearly incorrect at high frequencies, even though 
the Kondo resonance appears to be well resolved. For s — 10, 
the Hubbard peaks are still not resolved, but for s = 100, 
they are well reproduced. Finally, for s = 1000, the spec- 
tral function is fully resolved, including the band-edge fea- 
tures. We note that for s = 1000 the Matsubara points are 
chosen up to the asymptotic 1/z tail of the Green's function, 
see Fig fTOl c). These results are very instructive, since they in- 
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dicate that the information about the temperature is contained 
in the raw NRG results, not in the choice of the Matsubara 
frequencies. At low physical temperatures it is thus perfectly 
safe to use high fictitious temperature in the Pade approximant 
calculation. 
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Figure 1 1 : (Color online) Spectral function of the SIAM for different 
number A*'^ of the 2-averaging discretization meshes. For small N^, 
the spectra exhibit a large number of spurious resonances which are 
eliminated at larger A'^, . 



It is rather surprising that the Pade approximant produces 
smooth spectral functions with relatively little artifacts, given 
that the output from the NRG consists of a set of delta peaks 
at excitations energies which tend to be clustered. We indeed 
find that the Pade fitting applied to the results from a single 
NRG run with a specific discretization mesh (i.e., = 1) 
produces meaningless results, see Fig. [TT] However, already 
at A^^ = 2 the results are tremendously improved: the Kondo 
resonance, the Hubbard peaks, and the band-edge features are 
all resolved, located at the proper positions and with roughly 
correct spectral widths, although there are still sizable spuri- 
ous spectral peaks. With finer z-averaging, at iVr =4 and 
Nz ~ 8, the results are almost converged with only some mi- 
nor artifacts. At ~ 32 and beyond, we find no further 
improvement. 

In Fig. [12] we study the convergence with respect to the 
NRG discretization parameter A. The continuum limit is re- 
stored for A 1, but practical calculations are only possi- 
ble for A > 1.5. We find that the artifacts are indeed re- 
duced for smaller A, in particular in the Hubbard peak re- 
gion and nead band edges. The form of the Kondo reso- 
nance is well described by all A in the range considered, but 
there appears to be a small yet systematic trend toward higher 
Kondo resonance height for smaller A. We have verified that 
the Kondo resonance is adequately described even for much 
higher A = 4, but the description of the Hubbard peaks be- 
comes increasingly poor (results not shown). 
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Figure 12: (Color online) Spectral function of the SIAM for different 
values of the discretization parameters A. 
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Figure 13: (Color online) Shape of the Kondo resonance for T <^ 
Tk '■ close-up on the resonance peak and fits to various analytic ex- 
pressions. ISR stands for the inverse-square-root (Frota) form. 



C. The Kondo resonance 

We now discuss the shape of the Kondo resonance, which 
is well known to deviate strongly from the Lorenzian form, 
see Fig. [13] The behavior near w = is parabolic, as ex- 
pected for a regular FL systemi*^. The quadratic fit, how- 
ever, is only valid asymptotically for |ci;| <C Tk- The fit to a 
Lorentzian function is valid in a somewhat wider energy inter- 
val, but it starts to deviate appreciably already at a; ^ Tk- The 
long tails are better approximated by an inverse-square-root 



21.106.107 



as 



function (also known as the Doniach-Sunjic form) 
expected from the oithogonality-catastrophe physics for the 
scattering phase-shift 7r/2. The expression due to H. Frota 



,106 




2[1 + [lo/TkY] 



1/2 



(54) 
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Using the local-moment approximation method it has been 
shown that the tails are asymptotically logarithmi c'"^''"^ . 

At finite temperatures, the Kondo resonance can again be 
described by the Frota form using temperature-dependent pa- 
rameters. We write 



A{co,T) 



h{T) / 1 + ^1 + [uo/Tk[T)] 



1/2 



\^ 2{1 + [uj/Tk{T)Y) 
For the temperature-dependent height we find 

h{T) =\l + (2'/" - l){rT/TKY 



(55) 



(56) 



with s = 0.29, p = 1.67, and r = 0.66, while for the width 
we obtain 



TK{T) = cTK{l + a{T/TKf 



(57) 



with c = 0.60, a = 1.97, and b = 1.33. This expression 
well describes the Kondo resonance in the parameter range 
\uj\ < IOTk and T < lOT^. Note that we do not fix p = 2, 
as would be expected for a FL system. The T ^ asymptotic 
behavior is thus strictly speaking incorrect. Nevertheless, this 
choice gives a better description on the crossover scale T > 
Tk, which is more important for fitting experimental data. 

We test how well the Friedel sum rule ttTA{uj = 0) = 1 
for the particle-hole symmetric case is fulfilled in the low- 
temperature limit. The results for 7rrA(a; = 0, T) are shown 
in Fig.[T4]on the logarithmic temperature scale. The sum rule 
at T — > is fulfilled within 0.0018 using the broadening ap- 
proach, and within 0.0019 using the Pade approximant. The 
largest differences between the two approaches are in the in- 
termediate temperature region where the height of the Kondo 
resonance is systematically underestimated by broadening. 

Finally, we discuss an integrated spectral quantity, the 
conductance through a nanoscopic device (quantum dot) 
described by the SIAM, as given by the Meir-Wingreen 
formula^i^- 



G(T) = Go 



+ 00 



did 



5/ 
duj 



ttTA{lu,T), (58) 



where Go = ^ is the conductance quantum, and / is the 
Fermi-Dirac function. This quantity can be computed directly 
from raw spectral data, or using the continuous spectral func- 
tion obtain either by broadening or by Pade method. We find 
that the results agree very well in all three cases. We fit them 
on the phenomenological function proposed by Goldhaber- 
Gordon et al.ii^: 



G(r) = Go [l + (2^/^ - 1)(T/TkF 



(59) 



Note that G{T = Tk) = Go/2, which defines the Kondo 
temperature used in this work. This definition was used by 
Hamann, Nagaoka, and Suhl; it is sometimes denoted as 
Tk,h- In addition, there is Wilson's thermodynamic defini- 
tion x(T ~ Tk.w) / {giJ-BY' = 0.07, or the definition com- 
monly used in the perturbation theory works which is denoted 
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Figure 14: (Color online) (a) Kondo resonance height as a function 
of temperature, calculated using broadening and Pade approxima- 
tion, (b) Spectra for a range of temperatures from T /Tk = 10~^ 
to T /Tk = 10 in 13 equally spaced temperatures in the logarithmic 
scale. 



as T^. They are related through T^°^ = ^AYl^TKyv and 
Tk,h = '^■'2Tk,w- As discussed above, in a FL system, 
one may choose to fix p = 2 in order to obtain the expected 
T <C Tk asymptotics. For broadened spectra we find 



Tk = 4.4 • 10" 



0.229, 



while for Pade spectra 

TK==4.5•10"^ 5 = 0.227. 



(60) 



(61) 



The fitting was performed in the interval from T /Tk = 10^^ 
to T /Tk = 10. These values agree with the standard result 
s = 0.23^. By relaxing the constraint p = 2 a better fit 
is obtained in the intermediate temperature region where the 
cross-over occurs. In this case we find for broadened spectra 



Tk = 4.5 • 10" 



while for Pade spectra 



0.277, p = 1.74, (62) 



rK = 4.6-10"^ 5 = 0.279, p= 1.72. (63) 
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Figure 15: (Color online) Conductance G{T)/Go of a quantum dot 
described by the single-impurity Anderson model and a fit to the 
phenomenological Goldhaber-Gordon et al. formula. 



model). For a band of finite width, the Korringa-Shiba relation 
takes a slightly different formii^: 



lim 



Imx(w) 



2n [Rc((^,;53 + 5band))o] 



(68) 



where 5band is the spin of the conduction-band electrons, thus 
the correlator on the right-hand-side of Eq. (|68) is the im- 
purity magnetization induced by a small magn etic field ap- 
plied to all electrons in the systemii2. (See Ref. Il 13l for a re- 
lated discussion regarding the thermodynamic susceptibility 
and the Clogston-Anderson compensation theorem.) Alterna- 
tively, assuming proportionality between the magnetization in 
the impurity and in the band, one may work with Eq. (|67| i us- 
ing an effective parameter p. For instance, for non-interacting 
resonant-level model with e = and T/D = 10^^, one finds 
p w 0.889 X 2tt, while for T/D = IQ-^, p w 0.987 x 2tt. 
The correction to p/(27r) is thus approximately proportional 
to T/D. 



D. Dynamic susceptibilities 



The spin and charge susceptibilities are defined as 



Xc{z) = ^{{n;n))„ (64) 



where 5*^ = (l/2)(n^ — n^) and n = ri^ + n^, n„ being the 
impurity occupancy operators. We use units such that gfis = 
1, where g is the g ratio and /is is the Bohr magneton, and we 
set the electron charge e = 1. The factor 1/4 in Xc is added 
for convenience to make the two susceptibilities equal in the 
non-interacting limit U = for the symmetric case. 

For [/ = 0, the susceptibilities can be calculated exactly in 
terms of pair propagators 



2^ 



G„{uj + uj')G^, {uj')duj', 



— OO 

oo 



(65) 



as 



x.H = ^nf,M, xcM 



(66) 



The Green's functions here are time ordered, not retarded. 
The analy tic expression in the wide-band limit are given in 
Ref. Il 1 ll For a finite flat band they need to be computed 
numerically. Usually we plot the imaginary parts denoted as 
x"(w) = Im[x(w)], i.e., there is no — I/tt factor as in the 
spectral function. 

The Korringa-Shiba relation is a statement about the spin 
and charge dynamical susceptibilities in the low frequency 
limit. It can be written as 



lim^^.p[Rex(0)]^ 



(67) 



where p is a constant exactly equal to 27r in the wide-band 
Hmitliiii^ {T,U < £1 in SIAM, or J < D in the Kondo 




Figure 16: (Color online) Dynamical charge susceptibility of the 
non-interacting resonant-level model computed with broadening and 
with the Pade continuation. We also compare analytic results for 
the wide-band limit (D — > oo, dashed line) and finite flat-band (full 
line). T/D = 10"''. 

We first test the calculation of susceptibilities on the t/ = 
model, for a flat band with T /D ^ 10~^. Spin and charge 
susceptibilities are the same in this case, Xs(<^) = Xc(w). In 
Fig. [16] we compare the NRG results for x" obtained with 
broadening and Pade, as well as the analytic results in the 
wide-band limit and for the actual flat band. We observe 
that the susceptibility function has a peak associated with the 
charge fluctuation scale of F and a sharp drop at the band edge 
D. The analytic result obtained in the wide-band limit ob- 
viously does not account for the drop, which also explains 
the difference between the wide-band-limit coefficient 2tt and 
the effective coefficient p = 0.987 x 27r. The two numerical 
results overlap to a high degree, except near the band edge 
where the Pade approximant better describes the sharp de- 
crease. We find that the NRG result for the slope Imx(i^)/w 
agrees within one per mil with the exact result, while Rex(O) 
obtained via the Kramers-Kronig transformation of x" has a 
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one percent error. We thus obtain p = 0.967 x 27r, which is 
two percent lower than the exact result. Of course, Rcx(O) 
may be determined more accurately in NRG by other means. 



100 




T/D 



Figure 17: (Color online) Dynamical spin (top) and charge (bottom) 
susceptibility of the single-impurity Anderson model obtained with 
Pade continuation for a range of temperatures. The Kondo tempera- 
ture is Tk = 2.5 W'-^D. 



In Fig. [TT] we plot the dynamical spin and charge suscep- 
tibilities for the interacting SIAM in the Kondo regime for a 
range of temperatures. The charge susceptibility has a domi- 
nant peak on the scale of U (the maximum is at cj w 0.64C/), 
and a slight change of slope at the band edge oj = D. The 
spin susceptibility has a dominant peak associated with spin 
fluctuations on the scale of Tk, a hint of a spectral feature at 
the frequency where the charge susceptibility has a maximum, 
and a drop at the band edge D. With increasing temperature 
(in the range T ~ Tk), the charge susceptibility peak dimin- 
ishes slightly, while the main peak in the spin susceptibility 
shifts to higher frequencies and decreases in magnitude. We 
observe that the Pade approximant approach has difficulties 
with producing correct slopes in the lj ^ limit [see, for ex- 
ample, the T/D ~ 3 10^"^ curve in Fig. [T7l a)1. In order to 
determine the true origin of the problem, we studied the raw 
binned data and plotted the cumulant function 



which is relatively smooth for large 7V-. If x"{^) truly be- 
haved linearly near oj ~ 0, the cumulant should be quadratic 
for small x. We confirmed that this holds to a good approx- 
imation for X somewhat lower than the lowest energy scale 
of the problem (F for the non-interacting case, Tk for the 
SIAM), but not for x lower than the temperature T, where we 
find different power laws and the slope extraction becomes 
unreliable. This indicates that the Korringa-Shiba relation can 
only be tested in the T <C F and T <C Tk limits, respectively, 
but not for high temperatures where there is too much uncer- 
tainty. It also indicates that the raw dynamic information in 
the FDM NRG approach for w < T is not rehable. While this 
does not appear to be an issue for spectral functions, where the 
Pade approach produces what seems to be a reliable fit to the 
available raw data, this is not the case for dynamical suscepti- 
bilities where we observe incorrect slopes at low frequencies. 



IV. RESULTS FOR A CORRELATED-ELECTRONS 
PROBLEM: HUBBARD MODEL WITHIN THE DMFT 

We apply the Pade approximant approach to determine the 
spectral function of the Hubbard model in the paramagnetic 
phase within the DMFT. The questions of main interest are: 
i) Can we extend the upper limit of the temperature range 
where the DMFT(NRG) method can be safely applied towards 
T ^ D or even beyond? ii) Can we obtain more detailed in- 
formation about the internal structure of the Hubbard bands? 
iii) Is the problem of causality violation at lower temperatures 
reduced? 

We use the Bethe lattice with the non-interacting Green's 
function 



Go{z) = - 



z / z 

--zsgn(Imz)yi--^ 



(70) 



(69) 



The DOS of the Bethe lattice shares some features with the 
3D cubic lattice DOS. For instance, at the band edges it has 
square root singularities. The calculations are performed with 
Nz = 8 and we take advantage of the Broyden mixing to 
improve the convergence. 

We first discuss the half-filled system, (n) = 1, which 
is particle-hole symmetric. The local spectral functions are 
shown in Fig. [18] for two different temperatures. For low and 
moderate interaction U, the behavior is qualitatively the same 
as in the non-self-consistent SIAM. The Kondo resonance at 
low frequencies is reinterpreted as the quasiparticle band and 
there are two Hubbard bands which now exhibit some internal 
structure, in particular some enhancement at the inner band 
edges, which has also been resolved previously in the broad- 
ening approach with very narrow kernelsi^^ and is known to 
be a real feature of the DO S""^'"^ ; see, for example, the 
U/D = 2.5 curve at T/D = 10"^ Fig.[T8] 

The spectra at higher temperature, T/D = 10^^, are shown 
in the lower panel of Fig. [18] For values of U approaching the 
T = Mott metal-insulator transition, the quasiparticle peak 
is strongly suppressed and the spectral distribution inside the 
Hubbard bands is modified (see U /D ~ 2 spectra). Broad- 
ening seems to underestimate the height of the quasiparticle 
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Figure 18: (Color online) Spectral functions of the particle-hole sym- 
metric Hubbard model at half-filling, (n) = 1, at temperatures (a) 
T/D = 0.01 and (b) T/D — 0.1 for a range of repulsion strengths 
U/ D computed using the dynamical mean-field theory (DMFT). We 
compare spectra computed using the Fade approach (solid lines) and 
with broadening (dashed lines). (Fade continuation is used in all 
steps of the DMFT calculation, not just to obtain the final spectrum.) 



peak in this regime. For larger U/D = 2.5, the quasiparti- 
cle peak is completely washed out and the Pade results sug- 
gest that the Hubbard bands develop a pronounced three peak 
structure. This is an artifact of the method and broadening 
results show no such structure. 

In Fig. [19] we explore more carefully the reliability of the 
Pade approach in determining the internal structure of the 
Hubbard bands. In (a), we show a low temperature result 
where the Pade approach at A = 2 suggests, in addition to 
the well defined resonance at the inner band edge, also some 
feature at the outer Hubbard band edge. By reducing the dis- 
cretization to A = 1.7, this feature is washed out, which in- 
dicates that it was an NRG artifact. This is a general recepe: 
real spectral features may be distinguished from artifacts by 
changing the NRG parameters (such as A, truncation cutoffs, 
etc.), whereby real features should remain robust, while arti- 
facts are very variable. This is well illustrated in (b), where 
the results confirm that the three peek structure observed for 
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Figure 19; (Color online) Spectral functions of the particle-hole sym- 
metric Hubbard model at half-filling for U/D = 2.5: dependence on 
the NRG discretization parameter A. Dashed line show A{—lj)D. 



U /D = 2.5 at higher temperature is indeed an artifact. 

We next consider a hole-dopped Mott insulator with (n) = 
0.8 at U /D = 4. This same parameter set has also been used 
in a recen t study of the transport properties of this model, see 
Ref. Il 16i For low temperatures we find very good agree- 
ment between broadening and Pade approaches for the low- 
frequency part of the spectra, while at high frequencies the 
Pade approach reveals some internal structure of the upper 
Hubbard band which can only be obtained in the broadening 
approach if the broadening parameter a is much reduced. We 
note, however, that the integrated difference between two se- 
quential Green's function in the DMFT loop does not go to 
zero in the Pade approach but rather oscillates around IQ-^. 
This lack of true convergence can be shown to originate pre- 
cisely in the inner structure of the upper Hubbard peak, which 
changes subtly from iteration to iteration. Such behavior is 
sometimes associated with physical instabilities of the system 
which are not allowed for in the DMFT Ansatz, but may also 
here be an artifact. In any case, we believe that the upper Hub- 
bard band does have internal structure and differs from simple 
non-interacting DOS. 

We finally address the problems in the DMFT(NRG) at low 
temperatures which manifest as the violation of the causal- 
ity in the calculated self-energy functions, as discussed in the 
introduction. The self-energy has different asymptotic behav- 
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Figure 20: (Color online) Spectral functions of the Hubbard model 
describing a doped Mott insulator for a range of temperatures. We 
show curves obtained by broadening (a — 0.2) and Pade approxi- 
mation. 



ior than the hybridization function or the Green's function, so 
special care must be taken when performing continuation with 
Pade. For large frequencies, it behaves as^" 



S](z) - U{n)/2 + 0{l/z). 



(71) 



Because Pade requires 0(1/ z) asymptotic behavior, we first 
subtract U{n)/2, perform the continuation and add the sub- 
tracted value back to the result. It was believed that since the 
self-energy is a ratio of two Green's function and the real parts 



need to be obtained by the Kramers-Kronig procedure, it was 
supposed that the (over)broadening of the spectral functions 
(imaginary parts) at high frequencies has an impact on all fre- 
quency scales in the real parts, thus spoiling the causality. This 
appears not to be the case, since in the Pade approximant ap- 
proach the problems begin to show at the same temperature 
scale as with broadening. The causality violation thus appears 
to be intrinsic to the self-energy trick in NRG. 

Interestingly, we did not observe similar problems with the 
self-energy in calculations for the SIAM, but only in self- 
consistent DMFT calculations. 



V. CONCLUSION 

We have analyzed a non-standard approach for calculating 
finite-temperature spectral functions using the full-density- 
matrix numerical renormalization group (EDM NRG). It con- 
sists of evaluating the Green's function on the imaginary fre- 
quency axis, either at the Matsubara frequencies = i{2n+ 
1)ttT or at some other conveniently chosen set of points. To 
obtain the spectral function on the real frequency axis, an an- 
alytic continuation using the Pade approximant is performed. 
The technique works well in conjunction with sufficient z- 
averaging, which partially removes the discretization artifacts 
of the NRG. Compared with the broadening approach, the 
technique does not have any arbitrariness in the choice of the 
broadening kernel and the only free parameter in most cir- 
cumstances is the number of Matsubara points Nm', the low- 
energy part of the spectral function converges quickly with 
increasing N,n- In addition, unlike with broadening, the spec- 
tral function does not have artifacts on the frequency scale of 
w ~ T, which is particularly important for calculating trans- 
port properties (conductance, thermopower) of various sys- 
tems described by impurity problems, both in the context of 
nanodevices and in the dynamical mean-field theory (DMFT). 
Some artifacts in spectral functions appear on high frequency 
scales and we have shown that their origin is in the raw NRG 
data, not in the analytic continuation procedure as such. A 
further advantage of the proposed technique is its capability 
to resolve narrow spectral features away from the Fermi level, 
unlike in the broadening approach where the kernel width is 
usually proportional to the frequency and may thus far exceed 
the width of the spectral feature under consideration. The re- 
duction of overbroadening problems is sufficiently good that 
the self-energy trick of calculating the interaction self-energy 
function as a ratio of two correlators is not necessary. We have 
also shown that the technique may be applied at relatively low 
temperatures by choosing instead of the Matsubara points a 
modified set of points on the imaginary axis which provides a 
good sampling of the spectral information. We have tested the 
method first on the single impurity Anderson model, where 
we have analyzed the properties of the Kondo resonance both 
at low and at intermediate (T ~ Tk) temperatures. Since 
the method does not suffer from overbroadening nor from a 
lack of smoothness at finite temperatures, we were able to ob- 
tain the results for the Kondo resonance with unprecedented 
reliability. We discussed various analytic expression for the 
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Kondo peak shape and proposed an expression which is valid 
in a wide temperature range. We have also tested the approach 
on the Hubbard model within the DMFT. We find that the 
Pade approximant approach has a twofold advantage over the 
broadening; i) It has less spectral-function artifacts on the im- 
portant scale of a; ~ T. ii) It can better resolve the inner 
structure of the Hubbard bands. 

The main drawback of the Pade approximant are the arti- 
facts which, however, can be systematically eliminated by re- 
ducing A and increasing N^- In addition, we find that while 
the Pade approach works well for spectral functions, it has 
more difficulties with reproducing correct slopes of dynami- 
cal susceptibilities at finite temperatures (which is an issue of 
NRG itself). 

The importance of this work is not only in proposing an 
improved technique for obtaining finite-temperature spectral 
functions using Pade approximants, but more generally in 



showing that broadening of delta peaks is not necessarily the 
only nor the best approach. A possible improvement of our 
approach would be to build in as constraints additional infor- 
mation about the spectral function, such as the spectral mo- 
ments which are known to very high accuracy, and obtain 
the Pade approximant using an optimization procedure. One 
could also consider different fit functions. Work along these 
lines is in progress. 
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